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Abstract 

Cross-validation (CV) is often used to select the regularization parameter in high dimen- 
sional problems. However, when applied to the sparse modeling method Lasso, CV leads to 
models that are unstable in high-dimensions, and consequently not suited for reliable inter- 
pretation. In this paper, wc propose a model-frcc criterion ESCV based on a new estimation 
stability (ES) metric and CV . Our proposed ESCV finds a smaller and locally E'S'-optimal 
model smaller than the CV choice so that the it fits the data and also enjoys estimation 
stability property. We demonstrate that ESCV is an effective alternative to CV at a similar 
easily parallelizablc computational cost. In particular, we compare the two approaches with 
respect to several performance measures when applied to the Lasso on both simulated and 
real data sets. For dependent predictors common in practice, our main finding is that, ESCV 
cuts down false positive rates often by a large margin, while sacrificing little of true positive 
rates. ESCV usually outperforms CV in terms of parameter estimation while giving similar 
performance as CV in terms of prediction. For the two real data sets from neurosciencc and 
cell biology, the models found by ESCV are less than half of the model sizes by CV . Judged 
based on subject knowledge, they are more plausible than those by CV as well. We also 
discuss some regularization parameter alignment issues that come up in both approaches. 

Keywords: Lasso, model selection, parameter estimation, prediction. 



1 Introduction 



1.1 Regularization Methods 

There is an ever increasing amount of data in all fields of science and engineering. Often, this 
data comes in high dimensions relative to the sample size, posing a new challenge to scientists, 
engineers, and decision makers. These problems, plagued by the curse of dimensionality, suffer 
from overfitting when classical methods are applied. Regularization methods are used to tackle 
this problem of overfitting head on, usually by imposing a penalty on the complexity of the 
solution or through early stopping. For example, in fitting the usual linear regression model. 



the Lasso (Tibshirani 1996) and ridge regression (Tikhonov 1943 Hoerl 1962) adds a Li and L2 



penalty on the coefficient estimates respectively to the usual least squares fit objective function. 
Regularization methods can also take the form of early stopping iterative algorithms like classical 



forward selection or L2-Boosting (Friedman 2001 Biihlmann and Yu 2003 Zhang and Yu 2012 



Zhang 2011). Common to these methods is that they provide a family of possible estimators 



instead of just one estimator, with the unregularized solution at one end of the spectrum. This 
family is indexed by a regularization parameter and is commonly referred to as the solution path. 
For the Lasso and ridge regression, this regularization parameter determines the extent of the 
respective penalties. For the iterative algorithms, this parameter corresponds to the number of 
steps they take. Despite the difference in nature, numerous works have shown these regularization 



methods, at least in the context of the linear model, are intrinsically related (Efron et al. 2004 



Zhao and Yu 2007; Meinshausen et al. 2007). In that light, we will not focus on the distinction 



between the different types of regularization parameters but instead simply use A as a catch-all 
representation for them. In the same vein, we focus on the Lasso in this chapter even though we 
believe the method we present will work in the general framework. 



1.2 Selecting the Regularization Parameter A 

Much work has been done to show that regularization methods yield desirable solutions in high 
dimensional problems. For example, the popular Lasso has been shown to be L2-consistent 



(Meinshausen and Yu 2009 Bickel et al. 2009) and model selection consistent (Meinshausen and 



Biihlmann|2006 Zhao and Yu|2006 Tropp||2006 Wainwright|2009 ) in the high dimensional setting 



2 



when respective conditions are met. These results guarantee the existence of the A needed, but 
offer httle guidance on how to find the desired A in practice. Indeed, data-driven regularization 
parameter selection with guaranteed theoretical performance turns out to be a particularly difficult 
problem. 

One can rely on traditional model selection criteria like Akaike's information criterion (AIC) 



(Akaike 1974) and Bayesian information criterion (BIC) (Schwarz 1978). They are easy to com- 
pute but their validity rely on model assumptions. Furthermore, they are derived from asymptotic 
results, so even when model assumptions are satisfied, they may not work well in the finite sample 
case. 



More commonly used today are model-free approaches like cross-validation (CV) (Allen|[r974 



Stone 


1974 


) and bootstrap methods ( 


Efron 


1979 


Zhang 


1993 


Shao 


1996) 



computationally feasible for increasingly large data sets with the rapid advancements in computing 
power, especially the parallel computing paradigm that is currently the platform for dealing with 
big data. These methods rely on data resampling to assess prediction error of candidate solutions 
and can be found in various statistics and machine learning literature ( Hastie et al.||200^ Breiman 



1995 1996 2001). In particular, it is the most popular approach for regularization methods to 
select A. Doing so often leads to estimators with good predictive performance when the sample 
size is not small. However, there are other performance metrics that are also of interest in statis- 
tics, among them parameter estimation and variable selection metrics, with important practical 
connections. Unsurprisingly, optimizing predictive performance does not necessarily translate to 
having success with respect to these other performance metrics. 



1.3 Estimation Stability 

Statistical estimation is often tied to the optimization of an empirical loss or a random function 
based on data. Take for example, when fitting a linear model for random variables X G MP, Y & M, 
one might want to minimize the predictive L2 loss, 

f{f3) = Ex,Y{Y-X'^)'. 
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However, since the underlying joint distribution of (X, Y) is unknown, we instead minimize the 
empirical loss 

1 " 



n 

i=l 



where {xi,yi) for i = 1,. . . ,n, are the observed samples of {X,Y). By minimizing / instead of 
/, we incur a random estimation error dependent on the sample we observed. In the classical 
ideal scenario, e.g. when the sample contain independent and identically distributed observations 
and the sample size n is large and p is small, this estimation error incurred is small. If we draw 
multiple samples from {X,Y), each resulting estimate will be close to that of minimizing /, and 
consequently close to each other. This closeness across different samples can be seen as a form of 
stability in the estimation procedure, and we call it estimation stability. 

When the differences across different samples are measured by the L2 error, the estimation 
stability is obviously related to variance. We opt to use the term "stability" rather than the 
more commonly used term "variablity" in statistics. This is to recognize the fact that stability 
is a concept broader than variance or variability and that it is used in other quantative fields 
such as numerical analysis, dynamical systems, and linear analysis ( Higham|l996[ [Salle||1976 Ellis 



1998). Stability is also not associated with a particular metric (unlike variance) and thus allows 



its consideration under different metrics. In a recent paper (Yu 2013 (to appear), we advocate an 
enhanced emphasis on stability in statistical inference, especially for large and high dimensional 
data for which instability of statistical methods is much more common than in the domain of 
classical statistics. 

It is clear that estimation stability is a necessary property for a reasonable estimation proce- 
dure: the solution is not meaningful if it varies considerably from sample to sample. The converse 
certainly cannot be true in general: an arbitrary constant estimate will not vary but is certainly 
meaningless. When we add a data faithfulness requirement through cross-validation, we are able 
to devise a model-free criterion based on estimation stability for the selection of the regulariza- 
tion parameter A. Specifically, our proposed new criterion of estimation stability cross validation 
(ESCV) combines a new metric of estimation stablity (ES) with CV. For a given Li norm r or 
a regularization parameter A, our new ES{t) metric is the reciprocal of a test statistic for testing 
the null hypothesis that the regression function is zero. The test statistic is an estimate of the 
regression function standardized by an approximate delete-d Jackknife standard error estimate 
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based on the same pseudo data sets as in CV, and both estimates are functions of r. The pro- 
posed ESCV criterion chooses a local minimum of ES{t) which is smaller than the selection of 
r by CV. It is worth noting that the computational cost of ESCV is similar to that of CV and 
that they are both well suited to parallel computation, the dominant computing platform for big 
data. 

We demonstrate that our criterion ESCV provides a viable alternative to CV and BIC. We 
compare the three approaches with respect to several performance metrics when applied to the 
Lasso on both simulated data sets with different predictor dependence set-ups and two real data 
sets. These performance metrics are L2 error for parameter estimation, prediction error, F-measure 
and model size for model selection performance. We find that our criterion compares favorably 
with CV and BIC where they are known to excel, and outperforms them in other scenarios over 
different performance criteria. In particular, ESCV obtains excellent model selection results that 
are substantially better than those from CV , both in simulations and our real data sets, where 
the results are validated by subject knowledge. When the predictors are correlated, which is often 
the case in practice, ESCV also often outperforms CV for parameter estimation while at same 
time provides prediction errors comparable to those of CV . 

We note that previous works based on stability of solutions have shown positive results in 



terms of model selection (Breiman 1996 Bach 2008 Meinshausen and Biihlmann 2010). The 
work here differs from them in three substantial ways. Firstly, we develop a different measure of 
stability ES that is closely related to estimation than model selection, even though our ESCV 
does have desirable model selection properties such as the best F-measures across all simulation 
set-ups in Section [3j Secondly, we restrict our attention to selecting the regularization parameter. 
Even though we evaluate our choice by the performance of the corresponding solution, our focus 
remains on determining the right amount of regularization. We do not introduce any further tuning 
parameters as in Meinshausen and Biihlmann (2010). Concurrent with and independent of our 



work, recent follow-up papers on ( Meinshausen and Buhlmann||2010 ) use model selection stability 
to select edges in graphical models (Liu et al. 2010 Haury et al. 2012) or modify stability model 
selection to improve its false discovery rate theoretical properties ( Shah and Samworth||2013 ). The 
former two papers introduce further tuning parameters and they recommend fixed values for them. 



Shah and Samworth (2013) employs the complementary half-sample data perturbation scheme. 
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ESCV can work on such a scheme, but doing so would depart from the usual implementation 



of CV for comparison purposes. Thirdly, as in Meinshausen and Biihlmann (2010), these three 
papers apply data perturbation schemes such as bootstrap and subsampling with hundreds or 
thousands runs of model fitting. On the contrary, the CV (and ESCV) data perturbation scheme 
often works well based on 5-10 runs of model fitting. 



2 Methodology 

2.1 Lasso and Pseudo Solutions 

Let X G iR"^^, Y G M"' be our data set. The Lasso generates a family of solutions, 

/3[A]=argmm{||F-X/3||2 + A||/3||i}. 

/3[\], as a function of A > is also known as the Lasso solution path for f3j (j = 1, . . . ,p). We 
want to select a solution from this solution path; that is, choose a A and take its corresponding 
solution in the solution path. As alluded to earlier, we would like to make this choice based on 
estimation stability and fit. 

Since the notion of estimation stability is tied to the sampling distribution of the data, it is 
unavoidable that we need multiple solution paths to make such an assessment. Of course, it is 
often costly and infeasible to obtain extra data in practice. Thankfully, this problem is not new, 
and there are well-established ways to get around it. The key is to exploit the existing data 
by employing data perturbation schemes, parlaying it into multiple data sets. Let (X*[A;], F*[A;]) 
represents our kth pseudo data set, derived from (X, Y). In our case, these are the cross-validation 
folds: we randomly partition the data into V groups and form V pseudo data sets by leaving out 



one group at a time. (See Section 2.6 for other data perturbation schemes.) We then get pseudo 
solutions, 



P[k; A] = argmm - X*[m\l + A||/3||i} 



ioTk = l...V. 
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2.2 Alignment 

For many regularization methods, there are multiple representations for the regularization param- 
eter A. In the case of the Lasso above, A refers to the Li penalty parameter. Other popular choices 
to index the solution path are the Li-norm of the coefficient estimate, and the Li-norm expressed 
as a fraction of the Li-norm of the unregularized solution. Each of these representations for the 
solution path has its own merits, and is equivalent to the others (when non-trivial) for any single 
solution path. 

However, care must be taken on how to most meaningfully align our solution paths, when we 
reference the same A across different (pseudo) solution paths. In particular, when n < p, the Li- 
norm of the unregularized solution corresponds to the saturated fit and can vary a lot depending 
on which data points were sampled. This makes Li-fraction a poor choice, as the same index may 
correspond to very different amounts of regularization. The effect is more pronounced when the 
features are more correlated. Figure [T] shows a histogram of the maximum Li-norms for 10,000 
bootstrap Lasso estimates of the base case Gaussian simulation (with n = 100, p = 150, a = 1, 



p = 0.5) in Section 3.1.1 There is considerable spread: in this case, the upper decile is more than 
20% more than the lower decile. 

To highlight the effect of alignment on estimation performance, we compared the performance 



of cross-validation with the three alignments for the low noise scenarios detailed in Section 3.1.1 
As shown in Table [l| aligning the solution paths with Li-fraction does comparatively worse than 
aligning with Li-norm or the penalty parameter. Notably, in the popular R package "lars" used 
in solving the Lasso efficiently, the included cross-validation code aligns with Li-fraction. 

For ESCV to be proposed later, we find that there is little difference in performance when 
aligning with either the penalty parameter, A or the Li-norm r. In this work, since the Li norm 
of the solution is more comparable than the regularization parameter across different pseudo data 
sets, we opt to align with r. Note that our main comparison target CV performs best with the r 
alignment, as shown in Table [T] To be clear, 

r = r(A) = ||/3(A)||i, 

a one-to-one function for < A < Aq, where Aq = inf |a | ||/3(A)||i = o|. 

Indexing by r also has the benefit of visualizing the solution growing with the index. /3(r) 
starts at the origin when r = and moves towards the unregularized solution as r increases. 
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Figure 1: Empirical bootstrap distribution of maximum Li-norms of Lasso estimates on a typical 
simulated data set: a base case Gaussian simulation with n = 100, p = 150, a = 1, p = 0.5 in 
Section 13.1.11 

2.3 Convergence of Pseudo Solutions 

Given p-dimensional pseudo solutions /3[/c; r] for k = 1, . . . ,V , we want to measure their differences 
or see how similar or stable they are. Computing their pair-wise L2 errors was a natural first 
thing that we tried. However, we found that these errors vary too wildly to be useful even after 
normalization by means when there is high dependence between the components in the vector and 
this happens often especially when p is large. Notice that the components of an estimate of (3 
are combined in a linear fashion through X 13 to achieve our primary goal of estimating the linear 
regression function. Therefore we propose to compute the estimates 

Y[k-T]=X^[k;T], 

and study their stability. 

To evaluate such stability, as mentioned earlier we need a measure for how far apart the 
estimates are at each r: stable pseudo solutions should give similar estimates. One possibility is 



Histogram of max L1-norm 
of bootstrap Lasso estimates 
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Cross- Validation Estimation Error 


p 


Regularization parameter 


Li-norm 


Li-fraction 





0.794 


0.787 


0.817 


0.2 


0.781 


0.773 


0.821 


0.5 


0.969 


0.956 


1.04 


0.9 


1.83 


1.80 


1.95 



Table 1: Effect of alignment on cross validation performance on the base case Gaussian simulation 



with n = 100, p = 150, a = 0.5 in Section [3.1.1[ The first column corresponds to the alignment 
based on A, the second based on r and the third based on the Li fraction. Cross- Validation 
performs worst when aligning with Li-fraction. 



to look at the average pairwise squared Euclidean distance between the V estimates: 



A{r) 



Y[j;t] 



12- 



It is not hard to see that this is proportional to the more familiar "sample variance" formulation, 

V 



Var(r[r]) = l5^||r[fc;r]-f[r]"2 



k=l 

where f[r] = ^Er=i^[^;^]- 

Figure [2] shows two examples of this sample variance metric. The left panel is particularly 

illuminating: the pseudo solutions diverge as they grow at first but converge somewhat before di- 
verging again. Here, convergence and divergence simply refer to the sample variance metric (which 
is really just the average pairwise distance) decreasing and increasing respectively. Heuristically, 
this behavior is exactly what one would expect if there is a "correct" amount of regularization. 
Different samples would take different paths towards the "correct" solution before moving away 
from one another due to overfitting. Hence, we might select the r corresponding to the minimum 
point after the first negative slope. That is, we want to choose r corresponding to the "dip". 

By doing this, we incorporate fit into our selection even though our criterion is based on 
stability. The convergence of the solution paths is key: not only does it suggest we are close to the 
truth, we are also gifted with estimation stability. Note that this helps us automatically exclude 
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Figure 2: Examples of the sample variance metric. The left panel shows an example where the 
metric exhibits a "dip" , representing the "convergence" of the pseudo solutions. The right panel 
shows an example with a much muted "dip". It is difficult to use the sample variance metric to 
select a solution on the right. 

r's where the solution paths trivially agree. We see this trivial effect in Figure [2| where the global 
minimum for the sample variance metric occurs where the solutions are close to zero. 

However, this convergence effect is not always clear. The "dip" is not always present as shown 
in the example on the right panel. There you can still see the drop in gradient, but it is not clear 
which r we should pick. Notice, however, that in a solution path, the norm of the solution varies 
with the amount of regularization (by definition in our case). Since larger solutions naturally varies 
more, using the sample variance metric skews the choice towards solutions with small norms. We 
need to bring in the concept of normalization to account for this effect. 

2.4 Hypothesis Testing and the Estimation Stabihty Metric 

In hypothesis testing, a test statistic based on data is computed and its corresponding p-value 
is calculated by matching the test statistic with its model-specific theoretical distribution. This 
test statistic often takes the form of a mean value over its estimated standard deviation, e.g. the 
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student's t-test. The desired outcome for the t-test, as is often the case regardless of the assumed 
model and p-value computation, is to have thle test-statistic away from 0. The heuristic there is 
clear: if the hypothesized effect is real, the size of the mean value should be large compared to its 
estimated standard deviation. 

In the same vein, our sample variance metric should be relative to the squared mean size of 
the corresponding solution. We define the estimation stability metric, 

ES(r) := ^(M, 

the normalized version of the sample variance metric. Figure |3] shows the corresponding ES 
metrics in dashed lines superimposed on the old sample variance metric. On the left, the "dip" 
from the sample variance metric is preserved by the ES metric. On the right, there is now a 
pronounced minimum we can select. 





II -norm 



11 -norm 



Figure 3: Examples of the sample variance metric (as in Figure [2]) and the corresponding ES 
metric. We see that the ES metric preserves the local minimum from the sample variance metric 
on the left panel, and introduces one on the right panel where there was no local minimum from 
the sample variance metric. 



The ES metric's reciprocal has exactly the form of a test-statistic. We can view the ES 
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selection of r as a set of hypothesis tests. For each r, we are testing if the fit (^^[t]) is statistically 
different from fitting the null model {E{Y) = 0), albeit without a specified theoretical distribution. 
Our ES criterion of choosing the r corresponding to the convergence of pseudo solutions, is exactly 
choosing Y[t] with locally minimal normalized variance. This in turn, is exactly choosing the 
solution whose ES metric has the largest reciprocal, or in our analogy, the most statistically 
significant solution along the path. 



2.5 ESCV: Incorporating Cross- Validation 

There is no guarantee that our ES metric would have only one local minimum. Note that unless 
the multiple solution paths match up perfectly, there will be a local minimum or multiple local 
minima. Hence, even in the case where Y bears no relation to X at all, an inadvertent minimum 
on the ES metric will falsely suggest the pseudo solutions are converging towards a meaningful 
solution. To prevent scenarios like this where ES fails, we incorporate cross-validation into our 
selection. We have already limited our choice of minimum ES to local minima. Here we further 
limit it to the local minimum of r that is smaller than the cross-validation choice. We call this 
improved criterion estimation stability with cross validation (ESCV). In Section[3]on experimental 
results, we use a grid-search algorithm to find such a local minimum of ES as commonly done 
for CV. Thus ESCV's computational cost is similar to that of CV and they are both easily 
paralleable. 



We are exploiting the fact that cross-validation overselects (Leng et al. 2006 Wasserman and 



Roeder|[2009| . When gives a meaningful local minimum, cross-validation will likely overselect. 



Hence, ESCV behaves like ES above. However, when Y bears no relation to X, or when the 
noise overwhelms the signal, cross-validation will likely choose the trivial solution correctly. In 
this case, ESCV will follow suit and pick up the trivial solution. 

Note that this has negligible additional computation cost, as we are essentially getting the 
cross-validation choice for free. The bulk of the computation lies in computing the multiple 
solution paths we already have. 
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2.6 Discussion on ESCV 



Our ES metric is based on assessing the stability of the fitted values Y[t] = X(3[t] instead of the 
estimates /3[t]. This seems counter- intertuitive since we are interested in a variety of performance 
measures, most of which are based on the quality of /3[r] itself. However, we note that these 
performance measures only make sense if the underlying /3 is identifiable. To that end, there is a 
large volume of work showing the Lasso is model selection consistent under regularity conditions 
including that the smallest non-zero true parameter value is not too small compared a rate decay- 



ing in n (Meinshausen and Biihlmann 2006 Tropp 2006 Zhao and Yu 2006 Wainwright 2009). 



In particular, it assures us the asymptotic recovery of the underlying true (3 under appropriate 
conditions. 

However, in the finite sample case, and especially when the features are highly correlated, 
different linear combinations of features (of a given sparsity) may give approximately equivalent 
fits. Under data perturbation, it is not surprising that the different solution paths choose different 
features. This makes any metric based on /3[r] statistically unstable since V is small. Note that 
this does not contradict the assessment of the eventual f3[r] picked since ESCV and CV, picking 
from the same solution path, would both suffer from any failure of the original Lasso. 

In ESCV, we have used cross-validation folds to compute our pseudo-solutions. There are of 
course many other ways to generate pseudo datasets. One related approach would be to apply 



bootstrap sampling (Bach 2008). Here, simply sample with replacement from the original data 



set to generate multiple data sets. These two approaches are obvious choices, and can be applied 
to any estimation procedure (even those without an optimization formulation). A third choice, 
which applies only to penalized M-estimators such as the Lasso, is based on perturbations of the 



penalty (Meinshausen and Biihlmann 2010). Note that such perturbations of the penalty amount 



to perturbing (indirectly) the samples, but in a different way than bootstrapping. Finally, we can 
simply perturb the data directly by adding noise to X and/or Y. For example, we can add random 



Gaussian noise to the response (Breiman 1996). We find in the experimental results section that 



the choice of data perturbation scheme does not affect the results much. 

With high dimensional data, computation can be costly. In the case of the Lasso, the com- 



putation quickly gets expensive with larger data sets (Efron et al. 2004; Mairal and Yu 2012). 



Using the estimation stability metric to select the regularization parameter incurs only as much 
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computation as using cross validation. This is because the bulk of the computation in both cases 
rests in computing the solution paths of the V perturbed data sets. V in this case can be small 



as demonstrated in Section 3 This is in contrast to related work (Bach 2008 Meinshausen and 



Biihlmann 2010 ) which requires a much larger V. 



3 Experimental Results for Lasso 

In this section, we evaluate ESCV^s performance relative to the cross validation (CV) across a 
variety of data examples. In each problem, we fit a linear model using the Lasso. We focus our 
attention on the comparison with CV as it is the most popular criterion in practice. 

In all the data examples, we use the same grid-search algorithm to find a minimum of r for 
ESCV and that for CV. For our algorithm, we determine the domains for r, [0, Tmax] in each of 
our pseudo-solutions f3[l] r], where Tmax is the Li-norm of the saturated pseudo-solution. We limit 
our choice of r to the intersection of the domains. This potentially truncates some large values 
of r from the original lasso path, but we view this as inconsequential as we are often after sparse 
solutions corresponding to relatively small values of r. 

The fact that a r value is larger than the upper limit of the intersection implies that for at 
least one pseudo data set, the solution at this r value is not unique so this r value corresponds to 
instability and we would not want to consider it in our search for a stable solution. We evaluate 
each criterion in a equally spaced grid of 1000 values of r in the resulting domain. 

We start with simple sparse gaussian linear model simulations with our focus on the high 
dimensional data set up. We will vary the simulation parameters such as correlation strength 
within features and signal strength, as well as explore popular correlation structures of the design 
matrix, to cover a wide range of data scenarios in practice. We compare the solutions picked by 
ESCV and CV with regard to parameter estimation, prediction, and model selection performance 
measures such as F-measure and model size. We also include the BIC choice, but note that it 
performs poorly as expected in our high dimensional setting. 

We also explore the performance of our method on two real data sets from neuroscience and 
bioinformatics. We use a combination of objective predictive performance and subject knowledge 
on plausible models to illustrate the efficacy of ESCV over CV. In all cases, note that we are 
comparing different choices of r on the same solution path (from the original data). Furthermore, 
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we use the same data splits to make comparable results of CV and ESCV . We find that the 
number of cross validation folds V do not affect the relative performance between CV and ESCV . 
For the simulations and data below, we use commonly used \^ = 8 to obtain the pseudo solutions. 

3.1 Gaussian Simulation 

Let Xi G IBF for 2 = 1, . . . , n be independent identically distributed Gaussian variables with mean 
and covariance S. We have the usual linear model Yi = X^fS + ei, where /3 G MP is the unknown 
parameter, and G iR is independent Gaussian noise with standard deviation a. Pj are drawn 
from f/[|, 1] for j = 1, . . . , 10 and otherwise. The separation from zero is for model selection to 
make sense. This is a common assumption in theoretical work. 
The reported estimation and prediction errors are defined as 



respectively. For model selection, we use the F-measure which balances false positive and false 
negative rates of identifying non-zero coefficients of /3. The higher the F-measure the better. Each 
simulation is repeated 1000 times and the performance measures are aggregated across them. 

3.1.1 A Base Case 

Within the Gaussian linear model setup, there are many problem scenarios that favor one method 
over others. In particular, the following problem settings are known to affect the performance of 
the Lasso: correlation strength between features, strength of signal (size of coefficients) relative 
to the noise levels, dimension of the problem (p), and the correlation structure of the features. 
This is of course not an exhaustive list but is sufficient to cover a wide range of problems. As the 
strength of the correlation and signal are key to the behavior of the Lasso solution, we will include 
a full complement of these problem settings to illustrate when and why ESCV works well. 

We start with a base case scenario. Here, S has entries 1 down the diagonal and constant p 
on the off-diagonal. We vary p = 0,0.2,0.5,0.9 and a = 0.5, 1,2. We set n = 100 and p = 150 
to emulate the high dimensional data setting. Note that this implies that the columns of X are 
empirically correlated even when the features they represent are independent. 

As expected, CV does well in terms of prediction error (see Table[2]). However, observe that this 
does not necessarily translate to success in terms of other performance measures. With estimation 
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error, we find that once we leave the orthogonal case p = where estimation and prediction error 
are equivalent, ESCV has lower estimation error than CV despite having comparable prediction 
error. 

For model selection, we use the F-measure, the harmonic mean of the precision and recall 
rates, which are inversely proportional to false positive rate and false negative rate respectively. 
A high F-measure is achieved when both false positive and false negative rates are low. Recall 
that we are selecting solutions from the same solution path. The Lasso solution path corresponds 
roughly to a nested family of models in terms of features picked since features seldom gets dropped 
as the we relax the penalty term. Hence, having a low false negative rate (high recall) typically 
comes at the cost of a high false positive rate (low precision) . The F-measure balances these two 
objectives. 

By this measure, ESCV often outscores CV by a considerable margin. CV picks more true 
variables, but in the process picks up a disproportionately large number of noise variables. This is 



in line with theory that CV often overselects (Wasserman and Roeder 2009). ESCV cuts down 
the false positive rate, but not too much at the expense of the false negative rate. 

The results are summarized in Table [2] and the standard errors (SE) are given in Table [s] Note 
that the performance measures are highly correlated since for each simulation run, the selections 
by ESCV ^ CV and BIC are from the same solution path. Hence, the SEs for paired differences 
in performance measures are actually lower than the SEs for each of the values as reported in 
Table l 



3.1.2 Effect Of Ambient Dimension 

We repeat the simulations but this time for p = 50 and p = 500 to investigate the effect of the 
ambient dimension. Note that only the number of non-relevant features is changing; the number 
of non-zero coefficients remain at 10, the sample size n remains at 100. The comparison of ESCV 
and CV from the base case extends here: CV does well in prediction error, especially in the 
independent predictors case, but loses out to ESCV in the other scenarios with dependence more 
relevant to practice and in terms of parameter estimation and model selection metrics that are 
important for scientific applications. The results are summarized in Table [5] and |6| 

For the low dimensional case p = 50, we see that BIC performance is much improved compared 
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Estimation 


Prediction 


Model Selection 


Model Size 








error 






error 




_F-measurc 








p 


0" 


ESCV 


CV 


BIC 


ESCV 


CV 


BIC 


ESCV 


CV 


BIC 


ESCV 


CV 


BIC 





0.5 


0.470 


0.399 


0.632 


0.470 


0.399 


0.632 


0.673 


0.402 


0.225 


19.7 


39.7 


78.7 





1 


0.885 


0.787 


1.27 


0.885 


0.787 


1.27 


0.594 


0.413 


0.225 


21.7 


37.8 


77.4 





2 


1.53 


1.47 


2.29 


1.53 


1.47 


2.29 


0.434 


0.403 


0.207 


22.8 


29.8 


58.3 


0.2 


0.5 


0.386 


0.391 


0.687 


0.373 


0.360 


0.623 


0.535 


0.440 


0.223 


27.4 


35.4 


79.6 


0.2 


1 


0.757 


0.773 


1.36 


0.719 


0.711 


1.23 


0.522 


0.445 


0.223 


27.8 


34.5 


78.3 


0.2 


2 


1.37 


1.43 


2.56 


1.31 


1.33 


2.34 


0.475 


0.414 


0.206 


25.6 


31.6 


73.4 


0.5 


0.5 


0.465 


0.479 


0.871 


0.341 


0.348 


0.624 


0.493 


0.456 


0.222 


30.6 


33.8 


80.0 


0.5 


1 


0.927 


0.956 


1.72 


0.682 


0.695 


1.23 


0.485 


0.445 


0.214 


29.9 


33.6 


80.2 


0.5 


2 


1.62 


1.67 


3.18 


1.21 


1.23 


2.30 


0.427 


0.396 


0.19 


26.2 


29.5 


74.6 


0.9 


0.5 


1.02 


1.04 


1.92 


0.330 


0.339 


0.615 


0.466 


0.444 


0.211 


31.1 


33.2 


80.4 


0.9 


1 


1.75 


1.80 


3.55 


0.572 


0.587 


1.14 


0.396 


0.377 


0.183 


26.9 


28.9 


75.1 


0.9 


2 


2.46 


2.58 


5.50 


0.858 


0.861 


1.81 


0.274 


0.262 


0.158 


19.0 


21.3 


58.3 



Table 2: Performance of ESCV , CV and BIC in picking the regularization parameter for the 
Lasso for our base case design: constant correlation p, n = 100, p = 150. We see that ESCV 
performs best in parameter estimation (when different from prediction) and model selection, while 
doing comparably to CV in prediction. 



17 







Estimation 


Prediction 


Model Selection 






error SE 




error SE 




F-measure SE 


p 




ESCV 


CV 




ESCV 


CV 




ESCV 


CV 


BIC 


n 


u.o 


0.003 


0.002 


0.004 


0.003 


0.002 


0.004 


0.004 


0.003 


0.002 





1 


0.007 


0.005 


0.008 


0.007 


0.005 


0.008 


0.004 


0.003 


0.002 





2 


0.009 


0.008 


0.02 


0.009 


0.008 


0.02 


0.004 


0.003 


0.002 


0.2 


0.5 


0.002 


0.003 


0.005 


0.002 


0.002 


0.004 


0.002 


0.003 


0.002 


0.2 


1 


0.004 


0.005 


0.01 


0.004 


0.004 


0.009 


0.002 


0.003 


0.002 


0.2 


2 


0.007 


0.008 


0.02 


0.007 


0.008 


0.02 


0.003 


0.003 


0.002 


0.5 


0.5 


0.003 


0.003 


0.006 


0.002 


0.002 


0.004 


0.002 


0.003 


0.002 


0.5 


1 


0.006 


0.006 


0.01 


0.004 


0.004 


0.009 


0.002 


0.003 


0.002 


0.5 


2 


0.008 


0.009 


0.03 


0.007 


0.007 


0.02 


0.003 


0.003 


0.002 


0.9 


0.5 


0.006 


0.006 


0.01 


0.002 


0.002 


0.004 


0.002 


0.003 


0.002 


0.9 


1 


0.008 


0.01 


0.03 


0.003 


0.003 


0.01 


0.003 


0.003 


0.002 


0.9 


2 


0.01 


0.02 


0.08 


0.01 


0.005 


0.03 


0.003 


0.003 


0.002 



Table 3: Standard errors (SE) for performance numbers in Table [2] 
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to the p = 150 case. This is to be expected since BIC was developed to tackle problems in the 
classical regime of n > p. It is interesting to note that the performance of ESCV and BIC 
are very comparable in this case. Indeed, in this low dimensional case, they are very close to 
(empirically) optimal with respect to the F-measure. Of course, once we go to p = 500, BIC's 
performance falls off the cliff. ESCV continues to beat out CV by a wide margin. 

3.1.3 Other Correlation Structures 

The constant correlation structure can be seen as a simple one latent variable model. Here we 
introduce other correlation structures corresponding to more complex models and run the same 
simulations {n = 100, p = 150, and varying a and p). First, block correlation: all p features are 
randomly grouped into 10 blocks, and within each block, the features have correlation p while 
features from separate blocks are independent. Here, we let p = 0.3, 0.5, 0.9. Second, Toeplitz 
design: Ej^ = p''"-'', with p = 0.5, 0.9, 0.99. In this case, the ten true variables indices are randomly 
distributed among the p variables so that they are not all strongly correlated with each other. The 
results for the two designs are summarized in Tables [7] and [8] respectively. 

Despite the different correlation structures, the qualitative results from the prior section holds 
again in both variations. For prediction error, CV almost always outperforms ESCV, but ESCV^s 
predictive performance can be quite close to CV^s when p 7^ 0. For estimation error, ESCV gains 
on and eventually outperforms CV with increasing correlation levels. And for model selection, 
ESCV almost always has a higher F-measure than CV. Digging deeper. Table |9] shows the 
breakdown of the F-measure into the true positive and false positive rates. We can see that 
ESCV has much lower false positive rates while sacrificing relatively little on the true positive 
rates. Also unsurprising, BIC continues to do poorly in these p > n regimes. 

3.2 fMRI Data 

This data is from the Gallant Neuroscience Lab at University of California, Berkeley. In this 
experiment, a subject is shown a series of randomly selected natural images and the fMRI response 
from his primary visual cortex is recorded. The fMRI response is recorded at the voxel level, where 
each voxel corresponds to a tiny volume of the visual cortex. The task is to model each voxel's 
response to the n = 1500 images. The image features are approximately 10000 transformed Gabor 
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wavelet coefficients. We evaluate the prediction performance by looking at correlation scores 
against an untouched validation set of 120 images with 10-13 replicates. There are 1250 voxels in 
all. We ranked them according to their predictive performance under a different procedure from a 



previous study (Kay et al. 2008). Not all of them are informative, so we only look at the top 500. 

We find that while the prediction performance are nearly identical for ESCV and CV, ESCV 
selects much fewer features. The results are in Table |4} For the sake of brevity, they are averaged 
across groups of 100 voxels. For example, for the top 100 voxels, on average, the correlation 
scores are similar, but ESCV selects 30 features compared to CV's 70 features - a close to 60% 
reduction. That is, ESCV selects a much simpler and also more reliable model that predicts just 
as well as CV. Figure |4] shows how close the correlation scores are. 



Voxels 


Correlation Score 


Model Size 


ESCV 


CV 


ESCV 


CV 


1-100 


0.730 


0.735 


30.1 


70.2 


101-200 


0.653 


0.655 


27.0 


61.8 


201-300 


0.567 


0.566 


22.6 


49.6 


301-400 


0.455 


0.459 


16.7 


40.3 


401-500 


0.347 


0.347 


16.5 


33.6 



Table 4: Performance on fMRI data set. The numbers are averaged across the respective hundred 
voxels. ESCV cuts down the model size by more than half compared to CV, while largely 
preserving prediction accuracy. 



We note again that ESCV picks fewer features than CV by design (Section 2.5). That being 
said, the reduction is huge here: ESCV picks less than half the number of features as CV across 
the different voxels. Furthermore, this was with little or no loss in predictive performance. To 
understand the results better, we look at the individual voxels and examine the features selected. 
In almost all the cases, ESCV selects a subset of the features selected by CV. This is because 
they both select from the same Lasso solution path and features are rarely dropped after being 
added to the solution as we relax the regularization. 

Now, each feature corresponds to a Gabor wavelet characterized by its location, frequency, and 
orientation. We plot the features selected by both CV and ESCV as well as the extra features 
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Correlation Scores of ESCV vs CV for Top 500 Voxels 




0.0 0.2 0.4 0.6 0.8 1.0 
ESCV 



Figure 4: Scatterplot of predictive correlation scores of ESCV and CV for the top 500 voxels in 
the fMRI data set. We see that for almost all 500 voxels, the predictive performances are similar 
for ESCV and CV. 



selected by CV. The points in the plot represent the location and size of the Gabor wavelet 
selected. Figure [5] shows four randomly selected voxels. 

We can see quite clearly that the features selected by ESCV are clustered in one area whereas 
the features selected by CV but not ESCV are scattered across the image. Biologically, we expect 
each voxel to respond only to a particular area of the visual receptive field. This confirms that the 
extra features selected by CV are most likely not meaningful. Note that the location information 
of the Gabor wavelets were not used in fitting the model. 



3.3 Cytokine Data 



This data is from experiments performed by the Alliance for Cellular Signaling (AfCS), archived 
and made available at the Signaling Gateway, a comprehensive and free resource supported by 



the University of California, San Diego (UCSD). Pradervand et al. (2006) from the Bioinformatics 
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Figure 5: Feature selection by ESCV and CV on four randomly selected voxels. The "o"s 
represent features selected by both methods, while the "+"s represent features selected only by 
CV . The axes represent the pixel location of the images. The position and size of the points 
represents the wavelet location and wavelet scale respectively. Note that most of the extra features 
CV select are scattered and less biologically plausible. 



22 



and Data Coordination Laboratory at UCSD processed and analyzed this data in an attempt to 
identify signal pathways responsible for regulating cytokine release. There are 7 cytokines, 22 
signal pathway predictors. The signal pathways cannot be directly manipulated. Instead, ligands 
are stimulated to elicit responses from the signal pathway predictors and cytokines. For each 
cytokine, we have about 100 samples, each corresponding to average measured responses of the 
cytokine and signal pathways when a specific ligand pair is stimulated. 



In the original study (Pradervand et al. 2006), principal component regression (PCR) is used 
to fit the data to a linear model and select the significant signal pathways. The selection is done 
by thresholding the estimated coefficients via a pseudo-bootstrap method. They do this for each 
of the seven cytokines. That is, they solve seven linear regression problems, each with n ^ 100 
and p = 22, and apply thresholding to select the relevant signal pathways. These PCR results 
are then merged with other data and analysis to derive a final minimal model (MM). 

We run Lasso with ESCV and CV on the seven linear regression problems and compare our 
results with the results from PCR and MM. Fig [6] shows the feature selection results for the 
four methods. We regard MM as the benchmark for feature selection performance because it 
encompasses extra data and is not directly restricted by the linear model. 

We can see from Fig [6] that Lasso with CV does poorly. It selects the most features for every 
cytokine, often by a large margin. Lasso with ESCV on the other hand, selects the same or 
slightly larger number of features than MM. Moreover, with the exception of cytokine TNFa, 
ESCV always includes the features PCR selected which survived to the minimal model. In the 
case of TNFa, PCR barely selects (close to threshold) the one feature that ESCV missed. ESCV 
in general selects only about half the number of features PCR selects. There are far fewer false 
positives with respect to MM. At the same time, it rarely misses out any of the important features 
that PCR picked up. 



4 Conclusion 

Regularization methods are employed to deal with problems in the increasingly common high 
dimensional setting. However, the difficult problem of selecting the associated regularization pa- 
rameter for interpretation or parameter estimation, is not well studied. Our method ESCV is 
based on estimation stability but also takes into account model fit via CV. With a similar paral- 
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Figure 6: Feature selection results on cytokine data. The columns represent signal pathways 
predictors and each block of four rows correspond to a cytokine. The four rows within each 
block represent the selections of the four methods: the final minimal model (MM) and principal 
component regression [PGR) from the original study, and Lasso with ESCV and CV . The white 
squares corresponds to selected predictors. With only one exception, ESCV always selects the 
pathways that MM (which we regard as ground truth) does, while having much smaller models 
than CV. 
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lelizable computational cost as CV, we have demonstrated that ESCV is an effective alternative 
to the popular CV for choosing the regularization parameter for the Lasso. For the practical 
situation of dependent predictors, ESCV has an overall performance better than CV for param- 
eter estimation and significantly better for model selection. Their prediction performances are 
comparable, unless the predictors are independent. In particular, wc found much sparser models 
of less than half the size in both the real data sets from neuroscicnce and cell biology without 
sacrificing prediction accuracy, and these models are more plausible biologically based on subject 
knowledge. We beheve this result is not restricted to the Lasso but holds for other regularization 
methods as well. 

We also beheve that this method can also be readily extended to the classification problem 
through the generalized linear model, and leave this to future work. 
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Table 5: Performance of ESCV, CV and BIC in picking the regularization parameter for the 

Lasso for p = 50 with the base case Gaussian simnlation n = 100. constant correlation p. 
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4.85 
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1.56 
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26.9 
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Table 6: Performance of ESCV, CV and BIC in picking the regularization parameter for the 
Lasso for p — 500 with the base case Gaussian simulation n — 100, constant correlation p. 
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Table 7: Performance of ESCV, CV and BIC in picking the regularization parameter for the 
Lasso for the block correlation design, n = 100, p = 150. 
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Table 8: Performance of ESCV, CV and BIC in picking the regularization parameter for the 
Lasso for the Toeplitz correlation design, n — 100, p — 150. 
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Block design, p = 150 
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Table 9: Breakdown of the F-measure: the true positive and false positive rates of ESCV and 
CV for all the simulation scenarios. In all the cases above, there are 10 true variables and p — 10 
noise variables. 
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